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Abstract 


Topographic  Rossby  Waves  (TRWs)  have  been  identified  with  the  largest  variability  in  deep 
current  meter  records  along  the  continental  slope  in  the  Mid- Atlantic  Bight  (MAB).  Ray  tracing 
theory  is  applied  to  TRWs  using  the  real  bottom  topography  of  the  MAB  and  the  observed 
stratification.  The  dispersion  relation  for  TRWs  is  derived  and  various  wavenumber  limits  are 
discussed.  A  computational  method  for  tracing  the  waves  is  presented,  including  the  necessity  of 
smoothing  the  bathymetry.  In  the  examples  shown,  TRWs  with  periods  of  24-48  days  generally 
propagate  southwest  ward,  changing  their  wavelengths  from  400  to  100  kilometers  in  response 
to  the  change  in  bottom  slope.  TRW  paths  are  shown  that  connect  from  the  SYNOP  Central 
Array  near  68°W  to  the  SYNOP  Inlet  Array  near  Cape  Hatteras. 
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1  Introduction 


This  report  describes  the  method  of  ray  tracing  to  find  the  path  that  Topographic  Rossby  Waves 
follow  for  realistic  bottom  topography  and  stratification.  The  report  gives  some  background  on 
Topographic  Rossby  Waves  (TRWs)  and  explains  their  kinematics  and  dynamics.  It  then  documents 
a  set  of  programs  that  allow  for  the  development  of  TRW  ray  paths. 

2  Historical  Significance  of  TRWs  in  the  Mid- Atlantic  Bight 

TRWs  are  transverse,  quasi-geostrophic  waves  found  in  regions  of  sloping  topography.  They 
are  very  similar  to  ordinary  Rossby  Waves  in  that  the  sloping  topography  has  a  dynamically 
equivalent  effect  to  the  variation  in  Coriolis  parameter.  TRWs  thus  propagate  “westward”  along 
the  continental  slope,  traveling  with  the  shallow  water  on  their  right.  TRWs  have  been  detected 
by  moored  arrays  of  current  meters  at  numerous  locations  in  the  Mid-Alantic  Bight.  Their  strong 
horizontal  current  variance  overwhelms  the  signal  of  the  Deep  Western  Boundary  Current  and 
that  of  deep  expressions  of  Gulf  Stream  meanders.  By  measuring  the  horizontal  current  variance 
information  about  the  TRWs  is  obtained. 

The  source  of  the  TRWs  in  the  Mid-Atlantic  Bight  has  not  been  clearly  identified  yet.  One 
possibility  that  has  been  suggested  is  that  the  TRWs  are  created  by  large  Gulf  Stream  meanders 
and  ring  formation  [Hogg,  1981].  Although  there  appears  to  be  no  direct  coupling  between  TRWs 
and  the  above  Gulf  Stream  near  Cape  Hatteras  [Johns  and  Watts,  1985],  a  connection  between  the 
presence  of  TRWs  and  the  existence  of  troughs  and  ring  formation  downstream  has  been  reported 
by  Schultz  [1987].  Thus  it  appears  that  the  TRWs  may  be  remotely  forced  by  these  meanders  and 
rings  that  occur  farther  to  the  northeast.  This  remains  a  question  for  later  investigation. 


3  Dynamics  and  Kinematics  of  TRW  * 

To  derive  the  dispersion  relation  for  TRWs  one  solves  the  equations  of  motion  for  a  single  governing 
equation  in  terms  of  the  pressure  p.  The  equations  of  motion  consist  of  the  linearized  horizontal 
momentum  equations, 


du  _  _1_  dpt 

dt  V  po  dx 

dv  l  dpt 

dt  +  U  ~  po  dy 

the  incompressibility  condition, 

du  dv  dw  _ 
dx  +  dy  +  dz 

and  the  hydrostatic  equation  and  linearized  density  equation, 


(1) 

(2) 

(3) 


djf_ 

dt 


dp/  ,  A 


p0N2  w 
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(4) 

(5) 


These  are  then  expressed  quasi-geostrophically.  Because  TRWs  are  bottom  trapped  it  is  possible 
to  follow  Hogg  [1981]  and  assume  that  N  =  Nb{x,v)  the  Brunt- Vaisala.  frequency  near  the  ocean 
bottom.  The  governing  equation  becomes, 


d_ 

dt 


7, , ,  /o2  «v 


vv  + 


(6) 


Nb 2  dz 2 

Detailed  derivations  of  these  equations  are  given  in  Appendix  A.  The  relevant  boundary  con¬ 


ditions  are  no  normal  flow  through  either  the  top  or  the  (sloped)  bottom:  w  =  0  at  z  =  0  and 
w  =  u  •  —  Vh  at  z  —  - h(x ,  y).  These  conditions  result  in 


and 


ay  NB2wdh  dP‘  8M 

Wz  =  IT  ~  z=-h(z,y) 

Plane  wave  solutions  of  (6)  are  sought  of  the  form 


(8) 


p  =  A(z)  exp  [t(fcx  +  ly-  u>t)] 


(9) 


Substituting  (9)  into  (6)  and  using  the  two  boundary  conditions  (7)  and  (8)  it  is  possible  to  get 
the  dispersion  relation.  Specifically,  substituting  (9)  into  (6)  gives 

d*A 


dz2 


-  A 


£’(*.  +  !*  +  ») 
fo  w 


=  0 


(10) 


which  can  be  expressed  as 


where 


d2A 

dz2 


-  AM  =  0 


(11) 


Xi=(x),(t!  +  ,’+w>  <12> 

Here  A  is  the  bottom-trapping  coefficient,  not  the  wavelength.  The  solution  to  (11),  using  boundary 
condition  (7),  is 

A(z)  =  cosh  (Az)  (13) 


Substituting  (9)  into  the  bottom  boundary  condition  (8)  and  using  (13), 

Atanh(Aft)  =  hyk  -  hxl) 

JQU 


(14) 


Equations  (12)  and  (14)  constitute  the  dispersion  relation  for  TRWs.  It  is  convenient  to  non- 
dimensionalize  these  expressions  using  the  following  scales 


h  = 
A  = 


Hh * 

V 

H 


3 


T 


(--)  = 
Kdx'dy) 


0oj*NbH 

fo 

fo  (_d_  _d_\ 
NbH  \dz*’  dy* ) 


where  the  asterisks  denote  non-dimensional  variables.  The  factor  =  Rd  is  the  internal  Rossby 
radius  of  deformation.  The  non-dimensional  form  of  the  dispersion  relation  is 


A2  =  *2  +  /2  +  -  (15) 

u 

(l6) 

where  the  asterisks  have  been  dropped.  Appendix  A  shows  this  derivation  in  detail. 

Note  the  following  properties  of  TRWs: 

•  The  amplitude,  and  thus  the  energy,  decays  as  cosh(Az)  with  height  above  the  bottom;  i.e., 
the  wave  is  bottom  trapped. 

•  TRW  velocities  are  in  approximate  geostrophic  balance. 

•  TRWs  are  transverse  waves,  since  k  •  ug  =  kug  +  lvg  =  V  •  Ug  =  0. 

•  The  observed  periods  of  TRWs  in  the  Mid- Atlantic  Bight  range  from  8  to  60  days  [e.g.  Hogg 
1981;  Schultz,  1987]  for  wavelengths  of  40  to  300  km. 

The  dispersion  relationship  is  commonly  represented  either  <xs  two-dimensional  plots  of 

u>(fc)  or  as  “slowness”  curves,  which  are  traces  of  k,l  at  a  set  of  constant  frequencies  in  wavenumber 
space.  There  are  three  limits  of  (12)  and  (14)  that  will  be  examined: 

1)  Xh  -*  oo  and  correspondingly  tanh(Ah)  ~  1,  which  is  the  strongly  bottom-trapped  case. 

2)  Ah  — »  0  and  tanh(Ah)  ~  Ah,  which  is  the  uniform  “barotropic”  case. 

3)  0  =  0,  which  is  the  f-plane  assumption. 
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In  order  to  look  at  these  limits  it  is  necessary  to  choose  a  location  at  which  to  find  the  parameters 
that  will  be  used  in  the  dispersion  relation.  Pickart  and  Watts  [1990]  detected  TRWs  in  the  SYNOP 
Inlet  Array,  so  this  is  the  area  that  will  be  focused  on  when  discussing  the  limits.  Figure  1  shows 
the  location  of  the  SYNOP  Inlet  Array  and  sites  B1-B5  which  will  be  used  at  various  parts  of  the 
report. 

In  the  first  limit,  (14)  becomes 

A  =^-(hyk-hxl)  (17) 

JqU) 

which,  when  combined  with  (12),  gives  a  result  of  the  form 

Ak2  +  Bkl  +  Cl2  +  DI b  =  0  (18) 

where 

A  =  u2  -  NB2hy2 
B  =  2NB2hxhy 
C  =  -NB2h2 
D  = 

This  is  an  equation  for  either  a  rotated  hyperbola  or  a  rotated  ellipsoid,  depending  on  whether 
A  and  C  are  either  the  same  or  opposite  signs.  For  the  site  used  here  they  are  both  negative. 
Figure  2  shows  the  solution  to  (18)  for  a  wave  with  a  period  of  40  days  and  NB  =  .001s-1.  The 
bottom  slope  was  specified  as  hx  =  -0.018718  and  hy  =  0.0079937,  and  the  downslope  gradient  is 
arctan(^)  =  -23.125  deg.  These  parameters  were  chosen  because  they  correspond  to  site  B3  of 
the  SYNOP  Inlet  Array  (Figure  1)  where  TRWs  were  observed  by  Pickart  and  Watts  [1990].  This 
set  (a >,NB,hx,hy)  is  refered  to  as  the  “Test  B3”  values. 


5 


77 


75 


73 


71 


69 


67 


65 


Longitude  (W) 


Figure  1:  Location  of  SYNOP  Inlet  and  Central  Array.  The  circles  denote  the  sites  B1-B5.  The 
coastline  is  shown  by  a  dark  line.  The  other  dark  lines  are  four  TRW  traces.  The  trace  leaving  B2 
was  run  using  splines  fit  to  a  330  km  filter,  the  others  were  run  using  splines  fit  to  a  1 10  km  filter.  # 

The  B2,  B3,  B4,  and  B5  runs  are  14,  22,  33,  and  33  days  in  length,  respectively. 
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Figure  2:  The  bottom-trapped  dispersion  relation  for  a  TRW  shown  as  a  function  of  {k,l)  space. 
Units  are  m-1. 

In  the  second  limit,  tanh(Ah)  — *  (Ah),  (14)  becomes 

A  =  (t£)  -  '•*') 

When  this  is  combined  with  (12)  the  result  is 

(k  +  E )2  +  (l  +  Ff  =  G2 

where 

E  =  A_  IS. 

2u>  2 uh 


F= n± 

2u >h 


This  is  the  equation  of  a  circle.  The  solution  for  this  barotropic  case  at  site  B3,  assuming  the  Test 
B3  values,  is  shown  in  Figure  3. 


7 


In  the  third  limit,  0  =  0,  (12)  becomes 


which,  when  combined  with  (14)  does  not  yield  a  simple  analytic  solution.  Thus  this  limit  is 
solved  using  numerical  techniques  together  with  the  Test  B3  values.  Figure  4  shows  the  resulting 
dispersion  relation. 

In  reality,  TRWs  should  fall  somewhere  between  the  first  two  limits.  Figure  5  shows  a  numerical 
solution  for  the  dispersion  curve  using  the  Test  B3  values  but  without  specifying  any  limits.  A 
comparison  of  Figures  2  and  5  reveals  that  the  bottom-trapped  case,  where  tanh(Xh)  — *  1,  is  an 
excellent  approximation  for  the  full  solution  (at  least  for  waves  specified  by  the  Test  B3  parameters). 

Finally,  Figure  6  shows  the  dependence  of  the  full  dispersion  relation  on  frequency.  Numerical 
solutions  for  four  different  wave  periods  are  shown,  24,  32,  40,  and  48  days.  All  four  of  these 
wav  vectors  essentially  point  either  upslope  or  downslope;  the  best  alignments  were  found  with  the 
longest  period  waves.  As  seen  in  Figure  6,  for  wavelengths  greater  than  75  km  ( k  <  .0001),  k  and 
l  vary  only  slightly  for  these  moderate  changes  in  frequency. 

4  Ray  Tracing 

The  following  discussion  was  adapted  from  LeBlond  and  Mysak  [1978].  There  are  three  steps 
involved  in  ray  tracing.  First,  the  equations  describing  the  system  are  developed.  This  includes 
the  expressions  for  the  dispersion  relation,  the  group  velocity,  the  time  variation  of  the  frequency, 
and  the  dependence  of  the  frequency  on  the  environmental  parameters  (such  as  bottom  slope  and 
water  depth).  Second,  initial  values  for  position,  wavelength,  and  frequency  are  chosen,  and  the 
corresponding  group  velocity,  cff,  is  computed.  Third,  the  wave  packet  is  progressed  with  velocity 
c g  for  a  time  6t,  and  a  new  wavelength  and  frequency  are  found  by  integrating  the  equations  from 


Figure  3:  The  uniform,  barotropic  dispersion  relation  for  a  TRW  assuming  the  Test  B3  values. 
Units  are  m-1. 


Figure  4:  Numerical  Solution  for  the  0  =  0  dispersion  relation  in  which  the  Test  B3  values  are 
used.  Units  are  m~l . 
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Figure  5:  Numerical  solution  to  the  full  dispersion  relation.  Units  are  m-1. 
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Figure  6:  Numerical  solutions  for  periods  24,  32,  40,  and  48  days.  The  crosses  denote  a  period  of 
24  days,  the  square  is  a  32  day  period,  the  triangle  denotes  the  40  day  period,  and  the  diamond  is 
the  48  day  period.  Units  are  in  m-1 


step  1  over  the  time  interval  St.  The  second  and  third  steps  are  repeated  to  obtain  a  path. 

To  form  a  mathematical  expression  for  the  ray  tracing  consider  a  plane  wave,  with  amplitude 
A(x,t)  and  phase  5(x,  t), 

<£(x,t)  =  A(x,t)exp[iS(x,t)]  (19) 

It  is  assumed  that  the  amplitude  of  the  wave  varies  at  a  much  slower  rate  and  over  much  larger 
spatial  scales  than  the  phase.  The  local  wavenumber  k  is  defined  as 

k  =  VS 


and  the  frequency  u  as 


u  =  - 


ds 

dt 


Cross-differentiation  of  these  two  definitions  results  in  the  conservation  of  wave  crests  equation 


dk 

dt 


+  Vw  =  0 


(20) 


Here  k  is  the  directional  density  of  wave  crests  and  u  is  the  flux  of  wave  crests  past  a  fixed  point. 
By  considering  both  the  dynamics  of  the  wave  type  and  the  environmental  parameters,  it  is  possible 
to  relate  k  and  w  by  using  the  dispersion  relation 


w(x,t)  =  <r[k(x,t);7(x,0] 


(21) 


where  7  represents  the  environmental  parameters.  If  the  dispersion  relation  (21)  is  substituted  into 

the  conservation  of  wave  crests  (20),  the  result  in  tensor  notation  is 

dki  da  dk^  da  (hf  _ 
dt  dkj  dx{  +  dj  dxi 

Since  the  wavenumber  is  defined  as  the  gradient  of  the  phase,  its  curl  must  vanish;  thus  —  0 

and  so  Using  this  and  the  definition  of  the  group  velocity,  cg%  =  the  above  equation 

may  be  rewritten  in  vector  notation  as 


dk  =  9k 
dt  ~  dt 


do 

+  c*Vk=-^V7 


(22) 
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where  $  =  Jj  +  cg  •  V  is  a  derivative  moving  vrith  the  wavepacket.  An  expression  for  the  variation 
of  u  can  be  found  by  differentiating  the  dispersion  relation  (21)  with  respect  to  time  and  using  the 
conservation  of  wave  crests  (20)  to  eliminate  k.  The  result  is 


du >  _  do_  dki  do  d~{  _  do  du  do  d'f 
dt  dki  dt  +  &y  dt  ~  dk{  dx ,  d'l  dt 


Rewriting  in  vector  notation  similar  to  (22)  gives  the  variation  in  u  along  a  ray  path 


du  du  do  d'l 

dt  dt  W  d'l  dt 


(23) 


Thus  the  equations  which  describe  the  system  are 


w(x,<)  = 

dx 

dt 

dk 

dt 

du 

~di 


<r[k(x,t);7(x,f)] 

do 
C*  ~  dk 

for, 

d'l 

do  d'l 
d-f  dt 


(24) 

(25) 

(26) 
(27) 


5  Wave  Tracing  Programs 


This  section  describes  the  processing  steps  for  ray  tracing.  The  MATLAB  code  (or  FORTRAN, 
in  the  case  of  the  first  program  GET_ETOPO_II.FOR)  are  given  in  Appendix  C. 


5.1  Preparing  the  Topography 

One  of  the  most  important  features  of  this  method  is  that  real  bottom  topography  is  used.  In 
order  to  remove  all  small  scale  topographic  fluctuations  that  would  not  affect  the  waves  but  will 
influence  the  resulting  spline  fit,  the  topography  must  be  gridded  and  smoothed.  Figures  showing 
the  progressive  results  of  this  process  are  presented  in  Appendix  B. 


5.1.1  The  Raw  Bathymetry 


The  raw  bathymetry  can  be  obtained  from  the  NCAR  ET0P05  earth  elevation  data  set  which  is 
a  gridded  topography  array  on  a  grid.  The  data  set  is  gridded  in  equal  increments  of  latitude 
and  longitude  rather  than  kilometers  so  the  gridding  in  the  x  and  y  directions  is  slightly  different. 
Since  the  area  of  study  will  normally  encompass  only  a  small  fraction  of  this  data,  a  program  called 
GET_ETOPO_II.FOR  is  used  to  select  just  the  bathymetry  for  the  area  of  interest.  It  is  best  to 
choose  the  area  of  interest  generously,  so  that  boundary  effects  of  the  smoothing  filter  (discussed 
later)  will  be  minimized  in  the  actual  study  area. 

The  program  GET-ETOPO  JI.FOR  is  based  on  a  program  GET-ETOPO.FOR  (written  by 
John  Lillibridge  at  URI).  The  program  has  been  modified  to  remove  the  z-scaling  option  and  to  run 
on  a  micro  VAX  computer.  The  program  is  executed  by  the  command  file  GET_ETOPO.COM, 
which  must  be  modified  to  contain  the  appropriate  input  and  output  file  names. 

5.1.2  Remove  unwanted  features 

After  windowing  the  bathymetry,  it  is  necessary  to  remove  some  isolated  features.  In  this  for¬ 
mulation  only  waves  in  regions  where  WKBJ  approximations  are  valid  (i.e.  where  the  topography 
varies  on  a  much  longer  scale  than  the  wavelength).  Thus  features  that  violate  this  must  be  re¬ 
moved  or  the  subsequent  smoothing  process  will  spread  those  features  to  the  region  of  interset. 
For  example,  bathymetric  features  that  are  narrow  compared  to  a  wavelength  would  approximately 
generate  Taylor-columns  to  the  waves;  thus  it  is  better  to  eliminate  these  narrow  features.  The 
larger  features,  such  as  isolated  large-diameter  seamounts  or  chains  of  seamounts,  act  as  scatterers 
and  must  be  removed  since  reflection  is  not  included  in  the  formulation.  The  problem  features  are 
removed  by  reducing  their  vertical  extent.  For  example,  a  seamount  that  was  on  a  slope  of  average 
depth  -3500  meters  and  rose  to  -2800  meters  would  cause  problems  so  all  grid  points  in  the  region 


of  the  mount  that  had  a  value  of  greater  than  -3400  meters  would  be  set  equal  to  -3400  meters, 
effectively  shortening  the  seamount. 

The  MATLAB  routine  LANDSCAPE.M  is  used  to  eliminate  unwanted  features.  A  mouse 
is  required  to  run  the  program.  The  mouse  is  used  to  focus  in  on  the  feature  to  be  shortened.  A 
seamount  is  selected  and  eliminated  by  clicking  with  the  mouse,  first  on  the  lower  left  and  then  on 
the  upper  right  corner  of  the  seamount.  At  the  ‘landscape  level?’  query  give  the  desired  depth  of 
the  seamount  (e.g.  -3500).  The  program  will  continue  to  loop  through  these  steps  as  long  as  the 
user  types  ‘y’  to  the  ‘continue’  prompt.  The  bathymetric  contours  may  be  selectively  labelled.  This 
is  accomplished  by  clicking  with  the  mouse  on  the  desired  contour  at  the  ‘manual’  prompt.  After 
all  features  have  been  editted,  answer  ‘n’  to  the  ‘continue’  query.  Save  the  variable  ‘mowed_BATH’ 
to  a  disk  file;  this  is  the  new  bathymetry  data  set. 

5.1.3  Filtering 

After  removing  the  unwanted  large-scale  features,  it  is  necessary  to  filter  out  the  small-scale 
features  that  in  reality  do  not  affect  the  waves,  but  in  this  formulation  would  cause  considerable 
noise  in  the  calculation  of  the  environmental  parameters  7  and  V7.  The  TRWs  of  interest  have 
wavelengths  ranging  from  about  110  km  to  330  km.  Therefore  a  cutoff  wavenumber  of  was 
chosen  for  the  wavenumber  filter.  MATLAB  computes  the  filter  coefficients  from  the  ratio  of  the 
cutoff  wavenumber  to  the  Nyquist  wavenumber. 

The  program  for  filtering  the  topography  is  a  MATLAB  routine  called  FILTER2.M.  To  run 
the  program,  the  user  supplies  the  coefficients  for  the  spatial  gridding  and  the  ‘mowed-BATH’  data 
set.  The  smoothed  bathymetry,  Y,  should  be  saved  in  a  disk  file  for  later  use. 
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5.1.4  Compute  spline  representation 


The  final  step  in  preparing  the  topography  is  to  fit  a  series  of  B-splines  to  the  bottom  depths.  This 
is  done  so  that  the  gradient  of  the  depth,  V/i,  can  be  computed  for  any  given  location.  Because 
B-splines  are  linear  processes,  the  two  dimensional  grid  of  bathymetry  can  be  treated  one  dimension 
at  a  time.  The  B-splines  are  first  fit  to  lines  of  constant  latitude.  Then  another  set  of  B-splines 
is  fit  as  a  function  of  latitude  (i.e.  along  lines  of  constant  longitude),  creating  a  two  dimensional 
B-spline  expression  for  the  bathymetry. 

The  program  that  fits  the  B-splines  is  called  COMPUTE-SPLINES3.M.  The  code  is  com¬ 
mented  heavily.  Many  of  the  commands  are  concerned  with  reducing  “undulations”  at  the  edges 
of  the  domain.  These  are  reduced  by  augmenting  the  series  to  be  interpolated.  The  arguments 
of  the  program  are  the  latitude  and  longitude  of  the  first  point  of  the  filtered  bathymetry  data 
set.  It  should  be  noted  that  the  constant  ‘dlat’  must  be  specified  prior  to  executing  the  program. 
This  constant  specifies  the  grid  spacing  of  the  input  array  of  bathymetry.  (For  example,  for  a  grid 
spacing  of  half  a  degree,  set  dlat  =  0.5).  The  output  variables  ‘c,  cprime,  c2prime,  and  xknots’ 
should  be  saved  in  a  disk  file.  (Use  an  output  file  name  such  as  SPLINE_110  to  indicate  what 
filter  was  applied  to  the  input  bathymetry  data.) 

The  routine  TEST-SPLINES  uses  the  B-splines  variables  to  solve  for  the  depth,  bottom  slope, 
and  their  gradients  at  a  given  location.  It  is  executed  in  each  ray  tracing  step.  In  addition  to  the 
B-spline  coefficients,  the  latitude  and  longitude  of  the  point  must  be  supplied. 

5.2  Getting  the  initial  wavenumber 

As  noted  previously  in  step  2  of  the  Ray  Tracing  discussion,  it  is  necessary  to  obtain  initial  (Jt,,  /,) 
estimates  for  the  wavenumber  vector.  These  values  are  needed  in  order  to  compute  the  initial  group 
velocity  of  the  wave  packet.  The  MATLAB  routine  INITVAL.M  requires  the  user  to  supply  the 


position  of  the  starting  point  of  the  TRW,  0,  fo ,  Nb,  H,  and  Vh,  as  well  as  the  wavelength  and 
frequency  of  the  wave.  The  outputs  are  k  and  /  as  well  as  the  bottom  trapping  coefficient  A  and 
the  propagation  angle  9. 

Prior  to  running  INITVAL.M,  a  five  element  row  vector  called  ‘paras’  must  be  created.  The 
five  elements  of  this  vector  are  0,  dt,  fo ,  H ,  and  Nb-  The  second  element,  the  time  step  dt,  is  not 
used  by  this  program  so  it  is  not  crucial  to  supply  its  correct  value.  (However,  dt  will  be  used  in  a 
subsequent  program.)  The  third  and  fourth  elements  are  the  Coriolis  parameter  and  the  average 
water  depth  from  sea  surface  to  sea  floor.  The  final  element  is  Nb,  the  Brunt- Vaisala  frequency. 

The  routine  INITVAL.M  also  requires  the  depth,  h,  at  the  starting  point  of  the  TRW  as  well 
as  the  bottom  slopes  hx  and  hy.  To  find  the  value  of  h,  hx,  and  hy  at  this  position  the  user  must 
run  the  MATLAB  program  TEST-SPLINES.M  using  the  B-spline  coefficients  obtained  from  the 
COMPUTE-SPLINES3.M  program  together  with  latitude  and  longitude  of  that  point.  The 
output  of  the  TEST-SPLINES.M  program  should  be  placed  in  variable  called  ‘gradh’  for  use  by 
the  program  INITVAL.M  to  find  the  wavenumber. 

Four  additional  inputs  are  required  before  executing  the  INITVAL.M  program.  Set  ‘omega’ 
to  be  the  angular  frequency,  (3jr),  of  the  TRW  in  s-1  and  set  ‘LAM’  to  be  the  wavelength  in  km. 
The  longitude  and  latitude  of  the  starting  point  are  specified  as  ‘x’  and  ‘y.’  Once  all  of  the  variables 
necessary  to  run  INITVAL.M  are  specified,  the  program  may  be  executed.  The  initial  k  and  l 
values,  as  well  as  the  bottom-trapping  coefficient  A  and  the  propagation  angle  9  (measured  from 
east)  of  the  TRW  are  calculated  by  the  program. 

5.3  'fracing 

The  MATLAB  routine  called  TRWTRACE.M  traces  the  wave  path  in  the  manner  discussed  in 
the  section  entitled  Ray  Tracing.  The  routine  contains  a  self  check  to  see  if  the  error  is  becoming 
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too  large.  For  this  test,  a  new  frequency,  wj,  is  computed  using  the  new  wavenumbers.  This  is  then 
checked  against  the  initial  frequency  by  wcheck  =  .  Ideally  the  frequency  should  not  vary 

over  time.  However  due  to  computational  error  the  frequency  does  change  and  this  test  warns  if  the 
error  is  becoming  too  large.  Tests  (discussed  later)  show  that  for  the  40-day  period  waves,  errors 
in  wcheck  of  as  much  as  25%  will  not  significantly  alter  the  important  wave  parameters,  including 
the  path,  wavenumber,  group  speed,  or  the  bottom-trapping  coefficient. 

There  are  several  variables  that  need  to  be  assigned  prior  to  running  TRWTRACE.M.  First, 
the  time  step  increment  dt  (the  second  element  in  the  ‘paras’  array)  must  be  set.  This  program 
allows  the  wave  to  be  traced  forward  (positive  time  increment)  or  backward  (negative  time  incre¬ 
ment)  in  time.  The  time  increment  is  specified  in  hours.  Second,  put  the  longitude  and  latitude,  in 
that  order,  into  a  2-element  row  vector  called  ‘loc’  where  longitude  is  measured  in  degrees  east  (e.g. 
73W  is  -73).  Third  another  row  vector  called  ‘k’  should  containg  the  initial  ki  and  li  values.  Finally, 
the  user  must  specify  the  name  of  the  file  containing  the  B-spline  variables  (e.g.  SPLINE-110). 

This  is  an  approximate  guide  to  the  execution  time  of  the  TRWTRACE.M  program  on  a 
DEC  3100-76  computer.  The  major  time-consuming  step  is  in  the  calculation  of  Vh  which  requires 
the  computation  of  new  B-spline  variables  for  each  new  position.  The  execution  time  for  this  step 
is  lengthed  greatly  if  a  dense  spline  set  is  used.  For  example,  if  the  B-splines  have  been  fit  to  a  15 
by  25  grid  then  to  run  TRWTRACE.M  for  50  time  steps,  the  total  run  time  is  approximately  15 
minutes.  Denser  filter  grids  significantly  slow  down  the  program;  for  example  using  splines  fit  to  a 
85  by  145  grid,  a  50  time  step  run  would  take  1.5  hours. 

5.4  Examples 

Choosing  the  ray  tracing  parameters,  such  as  filter  size  and  time  step,  can  be  difficult.  One 
recommendation  is  to  conduct  test  runs,  specifying  different  choices  for  each  parameter  and  compare 
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the  results.  That  way  the  time  efficiency  can  be  weighed  against  the  desired  accuracy.  The  following 
examples  (Figures  6-10)  show  the  effects  of  various  choices  of  these  parameters  on  tracing  TRW 
ray  paths,  from  a  detection  site  in  the  SYNOP  Inlet  Array  to  a  potential  generation  region  in  the 
SYNOP  Central  Array.  (Thus  here  the  time  steps  will  be  negative). 

For  the  following  cases  a  consistent  set  of  parameters  was  used  and  only  one  parameter  was 
varied  at  a  time.  Unless  stated  otherwise,  the  bathymetry  wavenumber  filter  size  was  330  km,  the 
time  step  was  eight  hours,  the  period  was  40  days,  the  wavelength  was  130  km,  and  the  launch 
site  was  B3.  TRWs  of  this  wavelength  and  period  have  been  observed  in  the  Inlet  Array  (Pickart 
and  Watts,  1990).  Several  parameters  were  also  specified;  f3  =  1.8  X  10-11  m-1s-1  ,  /  =  9  x  10-5 
s-1,  and  H  =  4000  m.  As  mentioned  earlier,  Nb  was  set  equal  to  a  constant  value,  an  appropriate 
value  for  the  Gulf  Stream  region  being  Nb  =  1  X  10-3  s-1  (Pickart,  personal  communication). 

Figure  7,  shows  the  effect  of  varying  the  time  step  (dt)  on  the  wavepath.  Time  steps  of  one 
hour  and  eight  hours  were  tested.  The  paths  for  these  two  runs  are  almost  indistinguishable;  the 
only  difference  being  that  the  total  wave  distance  travelled  was  about  2%  shorter  for  the  dt  —  - 1 
hr  case  due  to  a  slower  group  speed.  Figure  8  compares  how  several  of  the  parameters  varied  over 
time  for  these  two  runs.  The  plots  of  the  total  wavenumber  (K),  the  average  bottom  slope  (gradh), 
and  the  bottom  trapping  coefficient  (A)  show  little  difference  between  the  two  runs.  However,  there 
is  a  very  large  difference  in  ‘wcheck’,  a  measure  of  the  percent  variation  of  the  actual  frequency 
in  the  run  compared  to  the  correct  frequency.  A  value  of  zero  for  wcheck  indicates  no  error  in 
computed  frequency.  In  Figure  8  it  is  obvious  that  the  frequency  errors  are  high  when  dt  =  —8 
hours;  however  this  appears  to  have  no  effect  on  the  predicted  wave  path  (Figure  6).  On  the  other 
hand  differences  in  the  actual  frequency  result  in  differences  in  the  group  speed  C„  for  the  two  test 
cases.  Figure  8  shows  that  when  dt  -  -1  hr  the  TRW  travels  about  2  —  3%  slower  than  when 
dt  =  -8  hrs. 
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Figure  7:  B3  launch  with  dt=— 1  hrs  and  dt=-8  hrs  for  410  time  steps  and  52  time  steps,  respec¬ 
tively.  Approximate  run  times  1.5  hours  and  15  minutes,  respectively.  The  cross  denotes  the  eight 
hour  time  step  and  the  square  is  the  one  hour  time  step  run.  The  dashed  lines  are  t.lje  330  km 
filtered  bathymetry. 
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The  effects  of  varying  the  degree  of  smoothing  of  the  bathymetry  were  also  tested.  For  example, 
for  the  area  between  the  SYNOP  Inlet  and  Central  Arrays,  the  110  km  filter  gives  a  more  realistic 
bottom  topography  than  the  330  km  filter  (see  Appendix  B).  However,  in  order  to  fit  B-splines  to 
the  110  km  filtered  bathymetry  without  aliasing,  the  filtered  data  can  only  be  subsampled  every 
|  degree.  This  leads  to  a  very  dense  rjline  set  and  considerably  slows  down  the  tracing  program. 
Using  a  330  km  filter  for  the  bathymetry  permits  subsampling  at  every  |  degree  without  aliasing, 
and  in  turn  leads  to  a  sparser  spline  set  and  a  shorter  run  time.  Figure  9  shows  the  predicted  wave 
paths  from  site  B3  for  three  different  filter  sizes,  110,  220,  and  330  km.  Despite  the  differences  in 
in  the  smoothed  bathymetries  the  end  positions  of  the  predicted  wave  paths  are  quite  close. 

The  influence  of  the  choice  of  launch  position  and  wave  period  was  also  tested.  Figure  10  shows 
launches  from  sites  B2,  B3,  and  B4  of  the  Inlet  Array,  which  are  about  50  km  from  each  other.  All 
other  parameters  have  been  held  constant.  The  three  paths  appear  to  parallel  one  another.  The 
biggest  difference  between  them  is  the  total  distance  traveled,  impling  that  the  group  speeds  are 
different.  The  slowest  speeds  were  obtained  for  the  wave  traced  from  site  B4,  where  |V/i|  is  the 
smallest. 

The  predicted  wave  path  also  depends  on  the  wave  periodicity.  Figure  11  shows  the  traces  of 
three  waves  launched  from  site  B3  with  periods  of  32,  40,  and  48  days.  In  this  case  the  waves  travel 
approximately  the  same  distance  (There  is  '  ..tie  variation  in  longitude  of  the  end  points).  On  the 
other  hand  there  is  a  difference  about  2°  latitude  between  the  trace  endpoints  for  waves  with  the 
32  day  periodicity  and  the  48  day  period  waves. 

These  results  can  be  summarized  as  follows.  Firstly,  an  accuracy  of  10  kilometers  in  the  ‘launch’ 
site  of  the  wave  gwes  consistent  path  predictions.  Additionally,  a  10%  error  in  wave  period  will 
also  give  consistent  path  predictions.  Second,  more  accurate  topography  is  obtained  when  a  110 
km  wavenumber  filter  is  applied.  Finally,  a  time  step  of  dt  =  -8  hours  is  adequate  for  most  cases; 
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Figure  9:  B3  launch  with  splines  on  110,  220,  and  330  km  filters  for  70  time  steps  of  dt=-8  hrs  on 
the  110  filter  and  52  time  steps  of  dt=-8  hrs  on  the  other  two  filters.  Approximate  run  times  1.5 
hrs.,  45  min.,  and  15  min.  respectively.  The  cross  denotes  the  330  km  filter  run,  the  triangle  is  the 
220  km  filter  run,  and  the  square  is  the  110  km  filter  run. 


Figure  10:  Launches  from  sites  B2,  B3,  and  B4,  each  for  52  time  steps  of  dt=-8  hrs.  Approximate 
run  time  15  minutes  for  all.  The  plus  denotes  the  B2  launch,  the  cross  is  the  B3  launch,  and  the 
circle  shows  the  B4  launch. 
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8  Appendix  A  The  Dispersion  relation  derivation 


This  is  the  derivation  of  the  dispersion  relation  (equations  (15)  and  (16))  that  was  shown  in  the 
Dynamics  and  Kinematics  of  TRWs  section.  Begin  by  using  the  same  five  equations,  the  two 
linearized  horizontal  motion  equations,  the  incompressibility  equation,  the  hydrostatic  equation, 
and  the  linearized  density  equation,  which  are  repeated  here. 


dt  po  dx 

dv  ,  _  1  dpi 

dt  +  U  po  dy 
du  dv  dw  _ 

dx  +  dy  +  dz 

dp'  poN2w 


Now,  taking  the  time  derivative  of  (31)  and  substituting  (32)  into  it  yields 


+  p0N2w  =  0 


Next,  take  ^  of  (29)  and  subtract  ^  of  (28)  to  give 

d2v  d2u  .(dv,  dv  \ 
dxdt  ~  dydt  +  f  [dx  +  dy)  +  ~  ° 

Substituting  (30)  for  the  quantity  in  paranthesis  and  moving  it  over  to  the  right  hand  side  gives 


d  (dv  du\  _  ,dw 
dt[dx  dy)+(3v  f  dz 


which  is  the  linear  vorticity  equation.  Let  the  quantity  in  paranthesis  be  denoted  as  Q.  Now,  by 
using  the  leading  order  geostrophic  balance  it  is  possible  to  get  relations  for  vg,  ug,  and  £g 


vg  = 


1  djf 
foPo  dx 
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(36) 


_ l_d£ 

Ua  foPo  dy 

■  s  — vV 

foPo 


Ca  =  7T-VV 


Here  V2  =  Substituting  these  into  the  linear  vorticity  (34)  gives  us  the  Quasi- Geostrophic 


Vorticity  (QGV)  Equation 


Now  take  of  (33)  to  get 


1  r  ^  /r-r2  A  ,  M  \  t  &w 

/opoWv  p)  +  (3dx^~fodz 


6  8  (  1  dj/\  dw 

dtdz  \poN2  dz )  dz 


which  is  substituted  into  the  right  hand  side  of  the  QGV  (38)  to  give 


1  .„2  A  ,  My  t  d  d  (  1  dj/\ 

•0p0  &(V  +  °dtdz  \poN2  dz) 


Multiplying  both  sides  by  p0  and  combining  the  time  derivatives  gives 


^.[VV  +  — [f —V  — l]  +  0—  -  o 

dt[  p±azl\NJ  dz"  P  dz 


Using  Hogg’s  [1981]  technique  approximate  the  varying  N  with  a  constant  Brunt- Vaisala  frequency, 
Nb,  at  the  bottom  to  obtain  a  workable  equation  of  one  variable. 

S  rr-rJ-/  .  f  fo  V  ^  .  M  _  „ 


IrW  +  ( AV  +  M.  =  n 

at^Vp+liVB/  dz^  +  ^dx  0 


This  is  the  same  expression  as  equation  (6).  Now  for  the  top  boundary  condition  at  z  =  0  the  rigid 
lid  approximation  is  used  to  force  w  =  0.  So  (33)  simplifies  to 


Assuming  no  normal  flow  at  the  bottom,  where  w  =  u  •  -  Vh 


-  _L 

~  fopo  [fly  dx  dx  dy\ 
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Substituting  into  (33)  gives 


ay  _  ni  r d£dh  _  d£m 

dtd 2  fo  [a®  dy  dy  dx\ 

where  z  —  -h(x,  y).  A  plane  wave  solution  of  the  form 


p'  =  A(z)  exp[i(fc*  +  ly-  wt)] 


is  sought.  Substituting  (42)  into  (39)  gives 


-  i2  +  ( A)  j  + 0ikP'  =  o 


From  this  point,  the  prime  will  be  dropped  from  pressure  p,  and  the  subscripts  will  be  dropped 
from  /  and  N.  Taking  the  time  derivative  yields 

f/V  la’A,  .  „ 


***  +  *-(£)  ;b?)  +  M?  =  0 


Divide  by  both  i  and  p,  since  they  are  non-zero,  leaving 


Dividing  by  u  and  regrouping  gives 


Multiplying  by  -(^)2A  produces 


S-Kt)'1*1*-*?']- 


If  the  quantity  inside  the  square  brackets  is  defined  as  A2  the  result  is 


—  -  A 2  A  -  0 
0*2  A  A  ~  0 


where 
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Using  (44)  and  the  top  boundary  condition  (40)  it  is  possible  to  solve  for  A(z) 

A(z)  =  cosh(Az) 


(45) 


Putting  this  into  (42)  and  then  substituting  into  the  bottom  boundary  condition  (41)  the  left  side 
of  (41)  becomes 


d2p 


dtdz 

and  the  right  side  of  (41)  becomes 


=  -Asinh(-Ah)ta;exp[*(A:x  +  ly  —  u>t)] 


m  i 

/  dh  dp 

dh  dp\ 

N 2 

/ 1 

Vdy  dx 

dx  dy) 

"  / 

)  =  ^j-i  cosh(-Ah)exp[t(fcx  +  ly-  u>t)]  k  - 


Setting  the  two  sides  equal  and  cancelling  the  common  i  and  exp[t(fcz  +  ly  —  u>f)] 

N 2 

-Asinh(— Xh)<j  =  —  cosh(—Xh)(hyk  —  hxl ) 

where  the  derivatives  of  h  are  denoted  as  hx  and  hv.  Dividing  through  by  cosh(-Ah)  gives 

N 2 

-A  tanh(-Ah)  =  ^j(hvk  -  hxI) 

Since  tanh(-x)  =  -  tanh(x),  this  can  be  written  as 

N 2 

A  tanh(Ah)  =  —:{hyk  -  hxl) 


(46) 


At  this  point,  the  equations  are  easier  to  work  with  if  they  are  non-dimensionalized.  Specify  the 
following  factors 


h 

A 


=  Hhm 


u  = 


y 

H 

f3u*NH 


(M)  = 
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Where  the  asterisks  denote  non-dimensional  variables.  Substituting  these  into  (44)  gives  the 
non-dimensional  form 


A*2  (N\2(km2f2  l*2f2  /3k*  f2  \ 

H2  ~  \f  )  \N2H2  +  N2H2  +  N2H2f3u • ) 


This  expression  can  be  simplified  by  multiplying  by  H2,  moving  ^5-  inside  the  parenthesis,  and 
cancelling  the  /3- terms  to  obtain 


A*2  =  Jfc*2  +  l*2  +  — 
w* 


(47) 


The  non-dimensional  form  of  (46)  is 


_  w1  /  si  (k  .*1 t  .  -LL' 

H  h(ff /  0U-NBNH  V*  NB~h’  NH. 


Cancelling  the  common  factors  gives 

...  r  (vt- - 

X  tanMA  A  )  -  — 

Finally,  introducing  the  radius  of  deformation,  Rd  =  yields 

...  .......  /  w-hsn 

X  Unh(A  A  )  = 


(48) 


Equations  (47)  and  (48)  together  comprise  the  non-dimensional  dispersion  relation  for  this  system. 
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9  Appendix  B  Example  Bottom  Topography  Preparation. 


Following  are  three  figures  that  show  the  progression  of  a  region  of  bottom  topography  as  it  is 
prepared  for  ray  tracing.  The  area  shown  is  65-77  W  and  34-41  N. 

The  first  is  the  original  bottom  topography,  as  taken  from  the  ET0P05  data  set.  Figure  12  still 
has  the  sea  mounts  and  all  of  the  true  topography  shown.  After  the  LANDSCAPE.M  program 
is  run  the  seamounts  that  were  in  the  upper  right  corner  are  no  longer  there,  creating  a  topography 
such  as  in  Figure  13.  Afte  e  filtering  has  been  done  the  result  is  Figure  14.  Here  a  110  km 
filter  has  been  used.  Figure  15  shows  the  220  km  filtered  topography.  Figure  16  shows  the  330  km 
filtered  topography. 


10  Appendix  C  Program  Codes 


Most  of  the  following  program  code  is  written  in  MATLAB.  A  good  reference  for  MATLAB 
is  Pro- Mat  lab™  for  VAX/VMS  Computers,  1991  by  MathWorks  Inc.  The  lone  exception,  the 


program  GET-ETOPO  JI.FOR  is  written  in  FORTRAN. 
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•OOTfflKJ 
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993 


(kf Inal , If Inal , laa, theta] -i 
a  [kf Inal , 1 final  ,1am,  theta] 
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